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Abstract 

We investigate the interplay between coherent effects characteristic of the propagation of linear 
waves, the non-linear effects due to interactions, and the quantum manifestations of classical chaos 
due to geometrical confinement, as they arise in the context of the transport of Bose-Einstein 
condensates. We specifically show that, extending standard methods for non-interacting systems, 
the body of the statistical distribution of intensities for scattering states solving the Gross-Pitaevskii 
equation is very well described by a local Gaussian ansatz with a position-dependent variance. We 
propose a semiclassical approach based on interfering classical paths to fix the single parameter 
describing the universal deviations from a global Gaussian distribution. Being tail effects, rare 
events like rogue waves characteristic of non- linear field equations do not affect our results. 



I. INTRODUCTION 



The progress on the experimental preparation and manipulation of interacting Bose- 
Einstein condensates has given a strong boost to the study of non-linear wave equations that 
account for the effect of interactions within the condensate in the framework of a mean-field 
approximation. Particularly promising cold-atom experiments in the context of transport 
physics include the realization of guided atom lasers [IHU , of arbitrarily shaped confinement 
potentials for cold atoms [5H7] , as well as of artificial gauge fields that break the time-reversal 
invariance for neutral atoms [HI E]- This makes it now feasible to experimentally explore the 
coherent transport of Bose-Einstein condensates through various mesoscopic structures that 
can possibly be modeled by billiard configurations. 

An interesting question that rises in this context is how the presence of the atom-atom 
interaction within the coherent matter waves affects interference effects well that are es- 
tablished for non-interacting systems. Indeed, previous theoretical studies have focused on 
the question how coherent backscattering in disordered potentials is modified by the pres- 
ence of the atom-atom interaction [TO]. These studies were recently complemented by our 
investigations on weak localization in the nonlinear transport through ballistic scattering 
geometries that exhibit chaotic dynamics [11] . While a semiclassical analysis of this non- 
linear scattering problem predicted a dephasing of the interference phenomenon that gives 
rise to coherent backscattering, signatures for weak antilocalization were obtained in the nu- 
merically computed reflection and transmission probabilities [11J. This effect was attributed 
to the specific role of non-universal short-path contributions, in particular to self-retraced 
paths the presence of which gives rise to a reduction of coherent backscattering as compared 
to the universal prediction. 

In the present work, we consider the same scenario as in Ref. [11], i.e., the quasi-stationary 
transport of bosonic matter waves through two-dimensional ballistic scattering geometries 
that exhibit chaotic classical dynamics. In contrast to Ref. [TT]^ however, we focus here 
not on transport observables such as the reflection and transmission probabilities through 
the billiard, but rather on the intensity distributions of stationary scattering states within 
the billiard. These intensity distributions are to be compared with the theoretical predic- 
tions that are obtained from the Random Wave Model (RWM) [121415] . which, in the linear 
case, represents probably the most powerful approach to describe the universal spatial cor- 
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relations of eigenstates arising from the classical chaotic behavior due to the presence of a 
spatial confinement. A most natural question that arises here is then to which extent the 
basic assumptions behind this model can also be used to describe possible universal spatial 
fluctuations in collective coherent matter waves that exhibit a weak atom-atom interaction. 
Within a mean-field semiclassical description, such matter waves are well described by the 
Gross-Pitaevskii equation [16] in which the presence of interaction is accounted for by means 
of a non-linear interaction potential. This equation is the starting point of our calculations, 
both on the numerical and on the analytical side. 

It is important to mention that rare effects due to the nonlinearity of the wave equation 
like rogue waves p2] or due to the presence of disorder, like branching [18j . will certainly 
affect the tails of the intensity distribution, and such effects are in principle outside the 
reach of our approach. Therefore, we will rstrict ourselves to the body of the distribution, 
where rare events need not to be considered. 

The paper is organized as follows. We describe in Section [TT] the scattering configuration 
under consideration as well as the main observable to be discussed in this work. In Section 



III we present a semiclassical theory of the intensity distribution in this nonlinear system, 
which is based on the Gaussian hypothesis as well as on the semiclassical theory of coherent 
backscattering in nonlinear systems. The predictions obtained by this semiclassical theory 



will be compared with the numerical results at the end of Section [HI] followed by a discussion 
in Section HVl 

II. STATIONARY SCATTERING STATES OF THE GROSS-PITAEVSKII EQUA- 
TION 

For our simulations, we use the inhomogeneous two-dimensional Gross-Pitaevskii equa- 
tion 

ih—^(r,t) = HV(r,t) + g(r) — \^(r,t)\ 2 ^(r,t) + S{r)e- lt » /h (1) 
where we have introduced the single particle Hamiltonian 

H= — [-ihV - qA(r)} 2 + V(r) (2) 
2m 

with the billiard potential V(r). This Gross-Pitaevskii equation contains a source term 

S(r) = S Xi(y)S(x - x L ) (3) 



which models the injection of atoms in a Bose-Einstein condensate acting as a reservoir with 
the chemical potential \x into the scattering system [12]. Xi(y) is a transverse eigenmode of 
the incident lead and So controls the current that is injected into the billiard. 

The non-linear potential term g(r)^\^>(r,t)\ 2 ^(r,t) describes atom-atom scattering 
events. Assuming that the degree of motion for the third spatial dimension is frozen out, 
e.g. by applying a harmonic confinement potential in this direction, we obtain the effective 
two-dimensional interaction strength as g(r) = 2y/2~7ra s / \JWj [mwi(r)] where a s is the s- 
wave scattering length of the atomic species under consideration and u± is the confinement 
frequency in the third spatial dimension. A spatial variation u± = Wj_(r) of this confinement 
will then naturally induce a corresponding variation in g = g(r). We shall, in the following, 
consider an effective interaction strength g(r) that is homogeneous within the billiard and 
vanishes in the attached leads. In a similar manner, we shall also assume that the artificial 
gauge field is given by A(r) = \Be^_ x r (with e_|_ the unit vector in the third spatial di- 
mension), with an effective "magnetic field" strength B that is constant within the billiard 
and vanishes in the leads. 

The billiard geometry considered in this work is shown in Fig. [1} It is characterized by 
the billiard area Q and the typical energy Eq of the incident matter-wave beam. Using these 
quantities, we can define a time scale h/ Eq, a length scale k^ 1 with Eq = fr 2 ^/ '(2m), and a 
scale i?o = 2nh/(qQ) (the flux quantum) for the magnetic field. All quantities in this work 
will be measured in these units. The area of the system is determined as ko^l 1 ^ 2 = 81.2. 
Two leads are attached to the billiard, which transforms it into an open scattering system. 
The width of the leads is given by W = 5Air/ko, which means that five channels are open 
in each of the leads. 

In order to calculate the stationary scattering states within this configuration, we insert 
the ansatz 

#(r,t) = m(r)e~ lt,l/h (A) 
into Eq. (jl|. This yields the self-consistent non-linear equation 

*(r) = S(r) (5) 

for the stationary scattering state. The amplitude of the source term is fixed such that 
in incident current of ji n = 1E /H is generated. Varying j in provides yet another way 
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FIG. 1: The shape of the billiard used in this work together with the density of a typical stationary 
scattering state. The hatched area in the left figure marks the region used for calculating the 
intensity distribution (adapted from Ref. [TT]). 

to effectively change the interaction strength g, as Eq. ^ is invariant under the scaling 
(gjin,®) ^ (gr]- 2 ,j in 7] 2 ,^r]) (for 77 e E). 

The non-linear scattering problem Eq. ^ is now solved using the methods described 
in Appendix [Aj We performed computations for 50 different values of the energy /i (all 
in the energy range 0.93 E . . . 1.18 E where five lead channels are open), for 25 different 
positions of the spherical obstacle in the centre of the billiard, and for the five different lead 
channels. The thereby obtained stationary scattering states \l/(r) are now used to determine 
the intensity distribution, i.e. the probability distribution of ^(r)! 2 , and its mean value. 
Only the points inside the marked region in Fig. [T] were used. Points in the vicinity of a 
boundary have to be avoided as explained in Sec. |III[ 

The left panel of Fig. [2] shows the probability distribution for obtaining a given real part 
of the scattering wavefunction (which is the same as for the imaginary part) within the 
marked region of the billiard in the linear (g = 0) and time-reversal invariant (B = 0) case. 
We find a very good agreement with a Gaussian distribution. Consequently, the intensity 
/ = l^p/d^/l 2 ) is distributed according to a Porter-Thomas law P(I) — e _/ , as is confirmed 
in the right panel of Fig. [2j There are tiny but systematic deviations from the Porter-Thomas 
law which slightly underestimates the actual intensity distribution near 1 = (as is also 
seen in the left panel of Fig. [2]) as well as for very large intensities / > 5, and overestimates 
it in between for 1 < / < 3. 

To highlight these deviations, we plot in Fig. [3] P(/)e 7 as a function of the intensity 
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FIG. 2: Left panel: numerical probability distribution (dots) of the real part of the wavefunction 
for B = and g = 0, which agrees very well with a Gaussian distribution (solid line). The 
same holds for the imaginary part. The right panel compares the numerically obtained intensity 
distributions P g {I) for g = and g = 0.05 (B = in both cases) with the Porter-Thomas 
distribution e~ J . Note the tiny but systematic deviations from the Porter-Thomas law, which 
are highlighted in Fig. [3] 

J, for various values of the nonlinearity g and the magnetic field strength B. A parabolic 
behaviour with a minimum at I — 2 is found. The prefactor of this parabolic scaling is 
reduced with increasing g. This appears natural as a weak repulsive interaction between the 
atoms is generally expected to give rise to a flattening of the density distribution, leading, in 
particular, to a significant reduction of intensity maxima, in order to minimize the interaction 
energy within the condensate (see also Ref. [20] for an analogous phenomenology in nonlinear 
optics in the presence of a defocusing nonlinearity). Indeed, similar findings were obtained 
for the quasi-stationary transport of Bose-Einstein condensates through two-dimensional 
disorder potentials [21J, for which is was found that the parabolic scaling of P(/)e 7 with 
the intensity / could even become inverted at stronger nonlinearities g. The dependence of 
the parabolic scaling with the magnetic field B, on the other hand, is related to coherent 
backscattering, for which we shall develop a semiclassical theory in the following section. 

III. THE SEMICLASSICAL APPROACH TO THE INTENSITY DISTRIBUTION 

In a first step, and following the now standard approach to describe the statistical proper- 
ties of eigenfunctions in non- interacting and classically chaotic billiard systems [12], we shall 
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FIG. 3: Deviation of the intensity distributions from the Poster- Thomas law for several values of 
the interaction strength g (left panel, with B = 0) and of the magnetic field B in units of Bq (right 
panel, for g = 0) 



make the fundamental assumption that scattering eigenstates of the non-linear Schrddinger 
equation share the same correlations as an ensemble of Gaussian Random Fields (see the 
left panel of Fig. [2J. This assumption leads to a Poster- Thomas distribution P(I) = e~ 7 
for the normalized intensity / = \^\ 2 / (|^| 2 ) (see Eq. ([7]) below) which, as discussed above, 
is supported by our numerical findings, as seen in the right panel of Fig. [2] The presence 
of a weak interaction does not change the excellent agreement of the numerical data with a 
Porter-Thomas profile. 

Knowing that the general features of the distribution of intensities for nonlinear waves 
are well described by a Porter-Thomas distribution, we now ask whether the deviations 
observed in Fig. [3] have also such universal character. Once again, the guiding principle 
will be linear case, where deviations from the body of the distribution are consistent with 
a Gaussian random field with a variance that smoothly depends on the local position. This 
consideration leads to an universal form of the deviations given by a Laguerre polynomial, 
which therefore depend only on a single parameter [13] . Fig. [3] shows how this property of 
the non-interacting case takes over perfectly when interactions are present. 
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The final step will be the explicit calculation of the coefficient in front of the polynomial 
corrections, and in particular its dependence on the strength of the interaction and of the 
magnetic field. Here we shall assume that a basic property of scattering states in the 
linear case, namely that their average intensity over energy and channels is related with the 
imaginary part of the full Green function, holds approximately in the presence of interactions 
as well. Assuming ergodicity within the billiard and utilizing the semiclassical approach 
presented in Ref. [11] , we obtain an explicit expression for the variation of the polynomial 
prefactor with the magnetic field strength for various values of the nonlinearity. 

A. The local Gaussian approach 

The calculation of the intensity distribution uses the values of |\l/(r)| 2 at many different 
positions, incoming channels, and energies. Thus, both an energy and a position average is 
involved. Motivated by the idea that for fixed position r, the average intensity over energy 
and channels (\^(r)\ 2 ^ E is itself a smooth function of r, the double averaging procedure is 
now split apart. 

We start therefore by assuming a position-dependent Gaussian distribution 

(6) 

for the real and the imaginary part of the wave function (\I/=\E' r + i^i) at a fixed point 
r, where (|\E'(r)| 2 ) £; denotes the energy and channel average of the intensity. For non- 
interacting systems with chaotic classical dynamics, such a local Gaussian distribution is 
a rigorous consequence of the Random Wave Model [22], and the possible universality of 
the deviations from the fully homogeneous case, i.e. from the case that {\^{r)^ E is inde- 
pendent of the position r, are encoded in (\^(r)\ 2 ^ E (see Ref. At a boundary the 
wavefunction vanishes and thus (^(r)) 2 ) vanishes there, too. Such boundary effects can 
also be incorporated in our approach. In this work, however, we shall restrict our study to 
the bulk, and therefore points in the vicinity of a boundary will be avoided. The distribution 
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FIG. 4: Comparison of the numerically ob- 



tained intensity distributions with Eq. (10). 



The unknown parameter f3 is determined by fit- 
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ting Eq. ( 10 ) (shown as dashed lines and marker 
symbols) to the numerical data (solid lines). 
The fitting is done in the range 1 < / < 7. 
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for the intensity p = \^\ is now calculated as 



Pr(p) 



Mr)\ 2 ) 
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Mr)] 2 ), 



(7) 



which is a local Porter-Thomas distribution for p. 

We now proceed by splitting (f)| 2 ) E into a position-dependent part and a position- 
independent part 

(8) 



1 



(Mr)\ 2 ) E = J [l + C(r)) 

by imposing the condition that the position average of C(r) is zero: (C(r)) = 0. Using 
(l) r = 1, we obtain relation ((\ty(r)\ 2 ) E ) = A^ 1 = (|\I/(r)| 2 ). Introducing the normalized 
intensity / = Ap, we can now rewrite the intensity distribution as 



Pr(I) 



1 



exp 



1 + C(r) 



-C(r)) 



exp 



-C(r)) 



-C(r)) 



e- 1 [C(r)]"L„(i) 



(9) 



n=0 



where in the last step we use the generating function of the Laguerre polynomials L n (I) [23J. 
Finally, we perform a position average to obtain, up to second order in J, the normalized 
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intensity distribution 



P{I)=(P r (T)) r = e- / [l + <C(r) 2 ) r L 2 (/)] 

1 + ( 1 - 21 + l -I 2 



(10) 

with f3 = (C(r) 2 ) r . In Fig. [4] we compare this formula with the numerically obtained 
intensity distributions. We see that the numerical data are very well described by a behaviour 



of the form (10), with (3 being the only free parameter. This supports our claim that for 



weak interactions, deviations of the intensity distribution are universal and depend only on 
a single parameter. 

B. Semiclassical calculation of j3 

The parameter can be numerically obtained by a fitting procedure and compared with 
a prediction based on the semiclassical approximation to the non-linear Green function 
G(r, r', E) defined through 



E-H-g(r)-\G(r,r',E)\ 
m 



G(r,r',E) =5{r-r') . 



(11) 



In order to understand the connection between the parameter (3 and the nonlinear Green 
function, we consider first the Green function Go for the linear system, 



[E-H] G (r,r',E) = 5(r-r'), 



(12) 



which admits a spectral decomposition in terms of the normalized scattering states ^E',a{ r ) 
at energy E' with incoming channel a, given by 



V ' ' 7 ^ J E-E'±i0+ 



(13) 



where + stands for an infinitesimal positive number. If we now consider the combination 
G+(r,r',E) - G (r,r',E) = / dE'^ E , >a {r)^ E ,, a {r')*8{E - E') (14) 

a 

we see that, up to numerical factors, 



(G+{r,r',E) - G (r,r',E)) E ex ^<* B)a (r)* s , a (r')*) B 



(15) 
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Therefore, the local variance ^|\I/(t*)| 2 ) can be calculated if we know the imaginary part of 
the Green function at r = r'. 

Although this construction depends on the fact that Gq is the Green function of a linear 
operator, our main assumption is that we can, for weak non-linearities, deform the linear 



scattering states into non linear objects such that a spectral decomposition of the form (13) 
for G holds, at least approximately. Following the same steps as for the linear case, we 
conclude that under such assumptions, the local variance for the interacting case is also 
related with the imaginary part of the nonlinear Green function. 

Although a closed expression for the non-linear Green function as a sum over classical 
paths is not known, it still satisfies a decomposition of the form 

G{r, r, E) = G zcm {r, r, E) + G long (r, r, E) (16) 

in terms of zero-length paths joining r with r' in zero time, and long paths hitting the 
boundaries several times. This decomposition carries over to the local variance which was 
defined in Eq. (8) as (\ty(r)\ 2 ^ E = [1 + C(r)\. The contribution from the zero-length 
paths produces then the uniform background l/A, while the long paths produce fluctuations 
around it to yield 

h 2 

C(r) = — rG lons (r, r, E) - G long *(r, r, E)} . (17) 
mi 

Finally, the average of C(r) 2 is computed within the diagonal approximation, where 
different paths are correlated only as long as they are related by time-reversal symmetry 
which is assumed to be weakly broken by the magnetic field. In perfect analogy with the 
derivation of the channel-resolved coherent backscattering probability that was calculated 
in Ref. [TTJ, we obtain 



P(B,g) 



mi 



G long (r,r)G lons *{r,r)) i 



r H r H l + {B /B w f x + ^.g (l + ^/b^ 



(18) 



All parameters in this formula are known. r# = mQ/h is the Heisenberg time, and the dwell 
time To as well as the characteristic scale B w for the magnetic field are determined by the 
classical dynamics of the system, as shown in Appendix |B| 



Figure 5 compares the semiclassical prediction ( 18 ) with the numerically determined value 



of f3. In the linear case g = the agreement is very good. In a similar manner as for the 
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FIG. 5: Comparison of the numerically determined values of f3 (see Fig. HI) in the left panel with 



the prediction from Eq. 18 in the central panel. Clear deviations are visible, but they have the 
form of a displacement (3 = /3 (g) that is roughly independent of the magnetic field and increases 
monotonically with g. In the right panel, this displacement /3 (g) is subtracted from the numerical 
data for (3 such that they match the prediction for B = ±0.5. 

channel-resolved retro- reflection amplitude [HI El], the parameter is enhanced for B = 
due to the constructive interference between trajectories that are backscattered from r to r 
and their time-reversed counterparts. Finite values of B introduce a dephasing between such 
trajectories, which leads to a suppression of the enhancement of the form ~ (1 + B 2 / B 2 ^)^ 1 . 



Eq. ( 18 ) predicts that the presence of a repulsive interaction gives rise to another dephasing 



mechanism for finite values of g, which, however, is slightly stronger for B = than for finite 
B and can therefore give rise, at finite but small values of g, to a local minimum of f3 (instead 
of a maximum) around B = (see the central panel of Fig. [5J. This minimum is found to 
be slightly more pronounced in the numerically determined values for (3 shown in the left 
panel of Fig. [5j As for the case of channel- resolved back-reflection [TTJ , this discrepancy can 
be attributed to non-universal short-path contributions, in particular to self-retraced paths 



whose contribution to (G long (r, r)G long *(r, r)) is doubly counted in Eq. (18) 



In addition to this minor discrepancy, we also find more significant deviations in the 
form of a global reduction of the numerical values for (3, which is independent of B and 
increases monotonically with g. Intuitively, this reduction could be explained by the general 
tendency of a defocusing nonlinearity to "smear out" the intensity distribution within the 
billiard, as was already mentioned above in the discussion of Fig. |3j Clearly, this tendency 
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would be independent of the presence of a magnetic field. A semiclassical evaluation of 
this effect, however, is beyond the scope of this work. It would, most probably, involve 
non-linear ladder-type diagrams that modify expectation values of higher moments of the 
local intensity as compared to the linear scattering problem. As we are, in this work, mainly 
interested in the dephasing behaviour of (5 as a function of the magnetic field, we subtract, 
in the right panel of Fig. [3j this global S-independent shift from the numerical data. Good 
agreement is then obtained with the semiclassical prediction. 

IV. CONCLUSION 

In this contribution we investigated, both numerically and analytically, the intensity 
distribution of non-linear scattering states. Our approach is based on a mean-field ap- 
proximation to the fully interacting problem of an atomic Bose-Einstein condensate, where 
interactions are incorporated by means of a non-linear term in the wave equation. Formally, 
we therefore expect that similar results hold for other kinds of non-linear wave equations, 
arising, e.g., in nonlinear optics. 

Our main finding is that not only the general features of the intensity distribution are 
universally reproduced by a standard Random Wave Model ansatz, but also that the small 
deviations from the body of the distribution can be understood in this framework by con- 
sidering local Gaussian statistics, in close analogy with the case of linear waves in classically 
chaotic geometries. We have finally shown that both the functional form of the deviations 
and their theoretical description by means of local modulations of the mean intensity are 
governed by a single numerical parameter. This parameter has an universal contribution 
originating from long ergodic paths which we were able to obtain in closed form by means 
of a semiclassical approach based on interfering classical trajectories. However, there is 
also a contribution that increases monotonically with the nonlinearity and is independent 
of the magnetic field, for which no theoretical approach is currently available. Once this 
latter contribution is identified and subtracted from the numerical data, we found very good 
agreement of the semiclassical approach with the exact numerical calculations. 
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FIG. 6: The black curve shows the transmis- 
sion obtained using the continuation method. 
% Results obtained through a time-dependent 

= i population of the billiard using Eq. @ are 

0.15 shown in gray. For large g, no dynamically sta- 
ble stationary solution exists. 

We would like to thank Josef Michl for valuable assistance in the evaluation of the semi- 
classical diagrams. This work was supported by the DFG Research Unit FOR760 "Scattering 
systems with Complex Dynamics" . 

Appendix A: Numerical computation of stationary scattering states 

In order to numerically solve the non-linear scattering problem, we discretize Eq. ^ 
using a second-order finite-difference approximation [2H]. This results in a two-dimensional 
irregular lattice whose lattice spacing is chosen such that we have approximatively 30 lattice 
points per wavelength. This ensures that the discretization error is negligible. The artificial 
gauge field A(r) is incorporated through a Peierls phase [26J. 

The interaction strength g(r) is assumed to be constant throughout the billiard but 
adiabatically ramped off inside the leads as explained in Ref. [27]. Therefore, the effects of 
the leads can be described, as in the linear case, by self-energies which provide the correct 
outgoing boundary conditions. This allows us to restrict the simulation to a finite spatial 
region. This procedure is analogous to the approach used in the recursive Green function 
method [281 EH]. 

The complex wavefunction \l/(r) is now represented by a 2A/"-dimensional real vector, 
with M the number of lattice points. Defining 

F : M 2Ar M 2Ar *(r) [^-fl]*(r) - g(r) — |^(r)| 2 ^(r) - S(r) , (Al) 

we have to seek for solutions of Ffif) = 0. This is done with Newton's iteration [30] 
^>k+i = — (T>F)~ l F(^k) which converges to a zero of F provided that the starting 
vector \l/o is suitably chosen. This choice is a non-trivial matter. Using g as an additional 
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FIG. 7: Probability distributions of path lengths L (left) and directed areas A (right) of classical 
trajectories inside the billiard that is shown in Fig. [TJ An exponential function is fitted (dots) onto 
both distributions after a short transient region. 



free parameter — i.e., g(r)=g go{f) with g e K and go{r) — 1 for r inside the billiard - 
we re-interpret F=F(^(r), g) as a function F : IR 2A/ " X R — > IR 2A ^. Neglecting critical points, 
the set F _1 (0) is a one-dimensional manifold which can be traced by a continuation method 
[301 13T] yielding the manifold as a parametric function s (^(s), g(s)) of the arclength s. 
An example of such a calculation is shown in Fig. [6j 

A prominent feature of non-linear wave equations is their potential multi-stability, i.e., the 
fact that they can support multiple solutions for a fixed value of g. In such a situation, the 
state that would be populated in an experiment depends on the history of the system. Here, 
we always use the the first solution found by the continuation method. This choice mimics 
the time-dependent population of the billiard that would be obtained from integrating Eq. ([T]) 
in the presence of an adiabatically slow increase of the source amplitude. 

Additional details of the numerical methods can be found in Refs. [HI 132] . 



Appendix B: Analysis of the classical dynamics 



The parameters td and B w in Eq. ( 18 ) can be determined using classical simulations. To 
this end, classical trajectories inside the billiard are calculated using a ray-tracing algorithm. 
The trajectories start in the left lead at a given longitudinal position x with a given total 
momentum p = > while the transverse coordinate y and the associated component 

p y of the momentum are randomly selected in a uniform manner. The simulation is continued 
until the trajectory leaves through one of the leads. 
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The time t 1 spent inside the cavity follows an exponential distribution 



Pit, 
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-t-t/TD 



(Bl) 



where td is the classical dwell time. Thus, an exponential fit (shown in Fig. [7]) of the nu- 
merically obtained path-length distribution yields the classical dwell time T£>. Its numerical 
value is, in our units, given by the average population ji n To ~ 241. 
A central limit ansatz results in the Gaussian distribution 



1 



cxp 



A 2 



2t 7 r] 



(B2) 



*/27r£^7 

for the directed areas A for paths of a given time t 7 . Here, r] is a geometry-dependent 
parameter that can be determined by evaluating the total distribution of A, 
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P(A) =t d P(t 1 ,A)e- t ^ /To dt 
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This is also an exponential distribution, and thus an exponential fit can be used to compute 
r) as shown in Fig. [7j The parameter B w is now finally determined as 
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We numerically find B w = 0.22 Bq in our units. Additional details are given in Ref. [TT 
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